数学建模
模型知识点
预测模型-线性回归-基本理论
Wang Haihua
🚅 🚋😜 🚑 🚔
线性回归¶
深入地学习线性回归之前，我们先来理清楚几个概念
- 确定性关系和相关关系：
- 参数之间的确定性关系是指某参数可以完全由其他参数所决定，例如正方形的面积与边长的关系，总价格与单价和数量的关系等。
- 相关关系是指一个或多个变量之间有一定的联系，但是这个联系不是完全决定性的。例如，人的身高与体重有一定的关系，一般来讲身高高的人体重相对大一些。但是它们之间不能用一个确定的表达式表示出来。这次变量之间的关系。我们称之为相关关系。又如环境因素与农作物的产量也有相关关系，因为在相同环境条件下农作物的产量也有区别，这也就是说农作物的产量是一个随机变量。
- 回归分析：回归分析就是研究相关关系的一种数学方法，是寻找不完全确定的变量间的数学关系式并进行统计推断的一种方法。它能帮助我们从一个变量取得的值去估计另一个变量的值。
- 线性回归：回归分析中最简单的是线性回归。
一元线性回归模型¶
单变量线性回归，又称简单线性回归（simple linear regression,SLR），是最简单但用途很广的回归模型。其回归式为：
$$ y=\alpha+\beta x $$为了从一组样本$\{\left(x_{i}, y_{i}\right)|i = 1,2,\cdots,n\}$之中估计最合适的参数$\alpha$和$\beta$，通常采用最小二乘法，其计算目标为最小化残差平方和：
$$ \min \sum_{i=1}^{n} \varepsilon_{i}^{2}=\min \sum_{i=1}^{n}\left(y_{i}-\alpha-\beta x_{i}\right)^{2} $$本课程我们暂不关心最小二乘法的推导过程，只需了解其原理及应用即可。感兴趣的同学可以参看这里
我们不加推导，直接给出最小二乘估计的结果
$$ \begin{aligned} \hat{\alpha}&=\bar{y}-\bar{x} \hat{\beta}\\ \hat{\beta}&=S_{x y} / S_{x x} \end{aligned} $$其中，
$$ \bar{y}=\frac{1}{n} \sum_{i} y_{i} $$$$ \bar{x}=\frac{1}{n} \sum_{i} x_{i} $$$$ S_{x x}=\sum_{i}\left(x_{i}-\bar{x}\right)^{2} $$$$ S_{y y}=\sum_{i}\left(y_{i}-\bar{y}\right)^{2} $$$$ S_{x y}=\sum_{i}\left(x_{i}-\bar{x}\right)\left(y_{i}-\bar{y}\right) $$现在有一组高中生的身高和腿长的数据，用一元线性回归研究他们之间的关系，然后对新的同学进行预测。先来看看数据的散点图
In [1]:
import numpy as np # 导入 numpy库，用于科学计算
import pandas as pd # 导入 pandas库 ，用于数据分析
import matplotlib.pyplot as plt # 导入 pandas库 ，用于数据可视化
%matplotlib inline
plt.style.use("ggplot") # 使用ggplot绘图风格
from sklearn.linear_model import LinearRegression # 导入线性回归工具函数 LinearRegression
x = np.array([143, 145, 146, 147, 149, 150, 153, 154, 155,
156, 157, 158, 159, 160, 162, 164]) # 输入x数据
x = x.reshape(16, 1) # 修改数据的格式，从行向量转换为列向量
y = np.array([88, 85, 88, 91, 92, 93, 93, 95, 96,
98, 97, 96, 98, 99, 100, 102]) # 输入y数据
plt.scatter(x, y) # 绘制散点图
plt.xlabel(r'$H$') # 添加xlabel
plt.ylabel(r'$L$') # 绘制ylabel
Out[1]:
Text(0, 0.5, '$L$')
很明显，升身高和体重明显符合线性关系。下面我们使用sklearn.linear_model中的LinearRegression方法进行线性回归
In [2]:
# 1. 初始化线性回归函数，命名为lrModel
lrModel = LinearRegression()
# 2. 使用lrModel对数据x,y进行拟合
lrModel.fit(x,y)
Out[2]:
LinearRegression()
In [3]:
# 3. 输出截距和斜率
print(lrModel.intercept_)
print(lrModel.coef_)
-16.07298072980727 [0.71935219]
In [4]:
# 4. 根据回归得到的系数，在散点图上绘制回归直线
LR_data = [lrModel.intercept_ + lrModel.coef_ * i for i in range(145,165)] # 求解回归直线上的一系列点
plt.plot(range(145, 165), LR_data) # 绘制回归直线
plt.scatter(x, y) # 绘制散点图
Out[4]:
<matplotlib.collections.PathCollection at 0x28149e849a0>
我们常常用相关系数$R^2$描述线性拟合的效果，使用方法如下
In [5]:
# 5. 计算R2
lrModel.score(x,y)
Out[5]:
0.928187845952738
执行代码可以看到，模型的评分为0.928,是非常不错的一个评分，我们就可以使用这个模型进行未知数据的预测了。（评分等于相关系数$R^2$用于表示拟合得到的模型能解释因变量变化的百分比，$R^2$越接近于1，表示回归模型拟合效果越好，通常认为$R^2 > 0.8$ 就是比较满意的回归结果了）
假定两位新同学的身高分别是170cm和163cm，下面使用本文的模型预测他们的身高。
In [6]:
# 6. 对新的数据进行预测
lrModel.predict([[170],[163]])
Out[6]:
array([106.21689217, 101.18142681])
多元线性回归模型¶
在回归分析中，如果有两个或两个以上的自变量，就称为多元回归。事实上，一种现象常常是与多个因素相联系的，由多个自变量的最优组合共同来预测或估计因变量，比只用一个自变量进行预测或估计更有效，更符合实际。
在实际经济问题中，一个变量往往受到多个变量的影响。例如，家庭消费支出，除了受家庭可支配收入的影响外，还受诸如家庭所有的财富、物价水平、金融机构存款利息等多种因素的影响，表现在线性回归模型中的解释变量有多个。这样的模型被称为多元线性回归模型（multivariable linear regression model）
$$ y=\beta_{0}+\beta_{1} x_{1}+\beta_{2} x_{2}+\cdots \beta_{m} x_{m} $$多元性回归模型的参数估计，同一元线性回归方程一样，也是在要求误差平方和最小的前提下，用最小二乘法求解参数。
$$ \min \sum_{i=1}^{n} \varepsilon_{i}^{2}=\min \sum_{i=1}^{n}\left(y_{i}-\left(\beta_{0}+\beta_{1} x_{1i}+\beta_{2} x_{2i}+\cdots \beta_{m} x_{mi}\right)\right)^{2} $$同样地，这里我们暂时不关心回归系数的具体的推导过程，只需要知道如何用Python操作即可。我们通过一个SO2浓度预测的案例进行介绍。
In [7]:
import pandas as pd # 导入 pandas库 ，用于数据分析
data = pd.read_csv('data/prediction_model/so2.csv') # 读入数据
print(data) # 显示数据
SO2(ppm) R G B S H 0 0 153 148 157 138 14 1 0 153 147 157 138 16 2 0 153 146 158 137 20 3 0 153 146 158 137 20 4 0 154 145 157 141 19 5 20 144 115 170 135 82 6 20 144 115 169 136 81 7 20 145 115 172 135 83 8 30 145 114 174 135 87 9 30 145 114 176 135 89 10 30 145 114 175 135 89 11 30 146 114 175 135 88 12 50 142 99 175 137 110 13 50 141 99 174 137 109 14 50 142 99 176 136 110 15 80 141 96 181 135 119 16 80 141 96 182 135 119 17 80 140 96 182 135 120 18 100 139 96 175 136 115 19 100 139 96 174 136 114 20 100 139 96 176 136 116 21 150 139 86 178 136 131 22 150 139 87 177 137 129 23 150 138 86 177 137 130 24 150 139 86 178 137 131
可以看到，这是一个五个变量的多元线性回归，我们先通过散点图来看一下各个变量单独与因变量的关系。
In [8]:
So2 = [0, 0, 0, 0, 0, 20, 20, 20, 30, 30, 30, 30, 50, 50, 50, 80, 80, 80, 100, 100, 100, 150, 150, 150, 150]
R = [153, 153, 153, 153, 154, 144, 144, 145, 145, 145, 145, 146, 142, 141, 142, 141, 141, 140, 139, 139, 139, 139, 139, 138, 139]
G = [148, 147, 146, 146, 145, 115, 115, 115, 114, 114, 114, 114, 99, 99, 99, 96, 96, 96, 96, 96, 96, 86, 87, 86, 86]
B = [157, 157, 158, 158, 157, 170, 169, 172, 174, 176, 175, 175, 175, 174, 176, 181, 182, 182, 175, 174, 176, 178, 177, 177, 178]
S = [138, 138, 137, 137, 141, 135, 136, 135, 135, 135, 135, 135, 137, 137, 136, 135, 135, 135, 136, 136, 136, 136, 137, 137, 137]
H = [14, 16, 20, 20, 19, 82, 81, 83, 87, 89, 89, 88, 110, 109, 110, 119, 119, 120, 115, 114, 116, 131, 129, 130, 131]
import seaborn as sns # 导入 seaborn库 ，用于数据可视化
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
colors = ['r', 'g', 'b', 'm', 'black'] # 指定每个绘图的颜色
xlist = ['R', 'G', 'B', 'S', 'H'] # 指定每个绘图的变量名
plt.figure(figsize=(20, 13)) # 指定图片大小
# 绘制R子图
plt.subplot(2, 3, 1)
sns.regplot(x=R, y=So2, color='r') # 绘制带有回归线和置信区间的线性回归图
title_text = 'R = {:.2f}'.format(np.corrcoef(R,So2)[1][0],)
plt.title(title_text) # 添加子图标题
# 绘制R子图
plt.subplot(2, 3, 2)
sns.regplot(x=G, y=So2, color='g') # 绘制带有回归线和置信区间的线性回归图
title_text = 'G = {:.2f}'.format(np.corrcoef(G,So2)[1][0],)
plt.title(title_text) # 添加子图标题
# 绘制B子图
plt.subplot(2, 3, 3)
sns.regplot(x=B, y=So2, color='b') # 绘制带有回归线和置信区间的线性回归图
title_text = 'B = {:.2f}'.format(np.corrcoef(B,So2)[1][0],)
plt.title(title_text) # 添加子图标题
# 绘制S子图
plt.subplot(2, 3, 4)
sns.regplot(x=S, y=So2, color='m') # 绘制带有回归线和置信区间的线性回归图
title_text = 'S = {:.2f}'.format(np.corrcoef(S,So2)[1][0],)
plt.title(title_text) # 添加子图标题
# 绘制B子图
plt.subplot(2, 3, 5)
sns.regplot(x=H, y=So2, color='black') # 绘制带有回归线和置信区间的线性回归图
title_text = 'H = {:.2f}'.format(np.corrcoef(H,So2)[1][0],)
plt.title(title_text) # 添加子图标题
Out[8]:
Text(0.5, 1.0, 'H = 0.83')
可以看到，除了第四个变量$S$以外，每一个变量基本上都和因变量有线性关系，因此我们可以删去$S$，进行多元线性回归。
In [9]:
from sklearn.linear_model import LinearRegression # 导入线性回归函数LinearRegression
lrModel = LinearRegression() # 初始化回归模型
lrModel.fit(np.array([R,G,B,H]).T,So2) # 输入需要回归的数据
print('截距为：', lrModel.intercept_) # 输出截距
print('系数为：', lrModel.coef_) # 输出系数
score = lrModel.score(np.array([R,G,B,H]).T, So2)
print('R2为：', score) # 输出相关系数R2
截距为： 2044.0215932049725 系数为： [ -1.38977229 -17.5022507 5.68341432 -9.34216398] R2为： 0.896131817285016
非线性拟合¶
多项式拟合（Polynomial regression）¶
尽管上述的线性模型已经能够求解许多实际问题，然而生活中的很多系统的增长趋势往往是非线性的，这时候我们就需要用到非线性拟合。
对于非线性回归问题而言，最简单也是最常见的方法就是多项式回归。多项式（Polynomial）是代数学中的基础概念，是由称为未知数的变量和称为系数的常量通过有限次加法、加减法、乘法以及自然数幂次的乘方运算得到的代数表达式。多项式是整式的一种。未知数只有一个的多项式称为一元多项式；例如$x^{2}-3 x+4$就是一个一元多项式。未知数不止一个的多项式称为多元多项式，例如$x^{3}-2 x y z^{2}+2 y z+1$就是一个三元多项式。
首先，我们通过一组示例数据来认识多项式回归，示例数据一共有 10 组，分别对应着横坐标和纵坐标。接下来，通过 Matplotlib 绘制数据，查看其变化趋势。
In [10]:
# 加载示例数据
x = [4, 8, 12, 25, 32, 43, 58, 63, 69, 79] # 输入x
y = [20, 33, 50, 56, 42, 31, 33, 46, 65, 75] # 输入y
plt.scatter(x, y) # 绘制散点图
Out[10]:
<matplotlib.collections.PathCollection at 0x2814ad66df0>
接下来，通过多项式来拟合上面的散点数据。一个标准的一元高阶多项式函数如下所示：
$$ y(x, w)=w_{0}+w_{1} x+w_{2} x^{2}+\ldots+w_{m} x^{m}=\sum_{j=0}^{m} w_{j} x^{j} $$其中，$m$表示多项式的阶数，$x^j$表示$x$的$j$次幂，$w$则代表该多项式的系数。当我们使用上面的多项式去拟合散点时，需要确定两个要素，分别是：
- 多项式系数$w$
- 多项式阶数$m$
这也是多项式的两个基本要素。
如果通过手动指定多项式阶数$m$的大小，那么就只需要确定多项式系数$w$的值是多少。例如，这里首先指定$m=2$,多项式就变成了：
$$ y(x, w)=w_{0}+w_{1} x+w_{2} x^{2}=\sum_{j=0}^{2} w_{j} x^{j} $$当我们确定$w$的值的大小时，就回到了前面线性回归中学习到的内容，具体来说，就是最小化残差平方和（最小二乘法）。
我们先尝试使用二次多项式拟合。
In [11]:
# 加载示例数据
x = [4, 8, 12, 25, 32, 43, 58, 63, 69, 79]
y = [20, 33, 50, 56, 42, 31, 33, 46, 65, 75]
from scipy.optimize import curve_fit # 导入非线性拟合函数curve_fit
# 定义需要拟合的函数形式，这里使用二次函数的一般式 y = ax^2 + bx + c
def f2(x, a, b, c):
return a * x**2 + b*x + c
plt.scatter(x, y) # 绘制散点图
popt, pcov = curve_fit(f2, x, y) # 执行非线性拟合
# popt数组中，三个值分别是待求参数a,b,c
y1 = [f2(i, popt[0], popt[1], popt[2]) for i in x] # 计算得到拟合曲线上的一系列点
plt.plot(x, y1, 'r') # 绘制拟合曲线
Out[11]:
[<matplotlib.lines.Line2D at 0x2814a7e8460>]
可以看到，其效果并不好，增大次数试一下？
In [12]:
# f3为三次多项式
def f3(x, a, b, c,d):
return a * x**3 + b*x**2 +c *x +d
# f4为四次多项式
def f4(x, a, b, c,d,e):
return a * x**4 + b*x**3 +c *x**2 +d*x + e
# f5为五次多项式
def f5(x, a, b, c,d,e,f):
return a * x**5 + b*x**4 +c *x**3 +d*x**2 + e*x +f
# 定义方差计算函数
def error(y1,y2):
a = np.array(y1)
b = np.array(y2)
return np.dot(a-b,a-b)
plt.figure(figsize = (10,10)) # 定义图片大小
plt.subplot(2,2,1) # 开始绘制第1张子图
plt.scatter(x, y) # 绘制（x,y）的散点图
popt, pcov = curve_fit(f2, x, y) # 执行2次多项式拟合
#popt数组中，三个值分别是待求参数a,b,c
y1 = [f2(i, popt[0],popt[1],popt[2]) for i in x] # 得到拟合曲线上的一系列点
plt.plot(x,y1,'r--') # 绘制拟合曲线
plt.title(str(error(y,y1))) # 计算方差，并作为图片的标题
plt.subplot(2,2,2) # 开始绘制第2张子图
plt.scatter(x, y) # 绘制（x,y）的散点图
popt, pcov = curve_fit(f3, x, y) # 执行3次多项式拟合
#popt数组中，三个值分别是待求参数a,b,c,d
y1 = [f3(i, popt[0],popt[1],popt[2],popt[3]) for i in x] # 得到拟合曲线上的一系列点
plt.plot(x,y1,'r--') # 绘制拟合曲线
plt.title(str(error(y,y1))) # 计算方差，并作为图片的标题
plt.subplot(2,2,3) # 开始绘制第3张子图
plt.scatter(x, y)
popt, pcov = curve_fit(f4, x, y)
y1 = [f4(i, popt[0],popt[1],popt[2],popt[3],popt[4]) for i in x]
plt.plot(x,y1,'r--')
plt.title(str(error(y,y1)))
plt.subplot(2,2,4) # 开始绘制第4张子图
plt.scatter(x, y)
popt, pcov = curve_fit(f5, x, y)
y1 = [f5(i, popt[0],popt[1],popt[2],popt[3],popt[4],popt[5]) for i in x]
plt.plot(x,y1,'r--')
plt.title(str(error(y,y1)))
Out[12]:
Text(0.5, 1.0, '38.052046270695655')
:::{admonition,dropdown,tip} 思考:一定是越高的阶次越好吗？
注意：兼顾拟合误差和预测误差。防止过拟合（overfit）和欠拟合（underfit）。
:::
一元其他非线性函数拟合¶
这里我们举一个指数拟合的例子与二次函数拟合进行对比。实际上，使用本课程讲到的方法，可以进行任意形式的函数拟合。
$$
f_1(x) = x^a + b
$$$$
f_2(x) = ax^2 + b
$$
In [13]:
# 定义f1
def f1(x, a, b):
return x**a + b
# 定义f2
def f2(x,a,b):
return a*x**2 + b
# 随机生成数据并绘制折线图
xdata = np.linspace(0, 4, 50)
y = f1(xdata, 2.5, 1.3)
ydata = y + 4 * np.random.normal(size=len(xdata))
plt.plot(xdata,ydata,'b-')
# 开始拟合
popt1, pcov1 = curve_fit(f1, xdata, ydata) # 使用f1进行拟合
popt2, pcov2 = curve_fit(f2, xdata, ydata) # 使用f2进行拟合
y1 = [f1(i, popt1[0],popt1[1]) for i in xdata]
y2 = [f2(i, popt2[0],popt2[1]) for i in xdata]
plt.plot(xdata,y1,'--',label = 'exp')
plt.plot(xdata,y2,'--',label = 'para')
plt.legend()
Out[13]:
<matplotlib.legend.Legend at 0x2814ab8ebe0>
